Figures

1 Import library

Code

# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)
library(gridExtra)

# table
library(gt)
library(gtable)

# map
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle and timezone
library(suncalc)
library(lutz)

# data wrangling
library(tidyr)
library(dplyr)
library(data.table)
library(purrr)
library(forcats)
library(lubridate)

2 Setting up custom function

2.1 windrose

Code

windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet_wrap = F,
           facet_grid = F,
           hour_sunset = 20,
           hour_sunrise = 7) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, hour_sunrise, hour_sunset), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin labels
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet_wrap == "origin") {
      data2$origin <- unique(data_to_plot$origin)
    } else if (facet_wrap == "trip") {
      data2$trip <- unique(data_to_plot$trip)
    }

    # add trip and origin if any
    if (facet_grid) {
      data2$origin <- unique(data_to_plot$origin)
      data2$trip <- unique(data_to_plot$trip)
    }

    # set up max value
    maxvalue <- 20

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "black",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet_wrap == "origin") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    } else if (facet_wrap == "trip") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(trip ~ .,
          strip.position = "right"
        )
    }

    # if facet
    if (facet_grid) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_grid2(trip ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

2.2 matlab_to_posix

Code

# function to convert matlab date into R posix
matlab_to_posix <- function(x, timez = "UTC") {
  days <- x - 719529 # 719529 = days from 1-1-0000 to 1-1-1970
  secs <- days * 86400 # 86400 seconds in a day

  return(as.POSIXct(
    strftime(
      as.POSIXct(secs,
        origin = "1970-1-1",
        tz = "UTC"
      ),
      format = "%Y-%m-%d %H:%M:%S",
      tz = "UTC",
      usetz = FALSE
    ),
    tz = timez
  ))
}

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code

# import the data
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartureLocation == "ANNU")

Code

# change few things
data_dive <- data_dive %>%
  mutate(
    # change divetype
    DiveType = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 100, 0, DiveType),
    DiveTypeName = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 100,
      "Transit",
      DiveTypeName
    )
  )

# add Post Molting vs Post Breeding trip
data_dive <- data_dive %>%
  group_by(id) %>%
  mutate(trip = if_else(between(month(first(date)), 1, 3), "Post Breeding", "Post Molting"))

4 Data Visualisation

4.1 Table 1

We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.

Code

compute_table_1_a <- function(x) {
  # if the sum of the column is below 100
  if (sum(x) < 100) {
    # it's definitely a proportion
    # so we weighted average these
    y <- weighted.mean(x, c(1068746, 2207361))
    y <- paste0(round(y * 100, 1), "%")
    return(y)
    # otherwise we just sum the column
  } else {
    y <- prettyNum(sum(x), big.mark = ",")
    return(y)
  }
}

Code

# shark and orca area = NA if no location data
data_dive <-
  data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
table_1_a <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T)
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives per area",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark",
    n_benthic_dives_in_orca_proportion = "Orca",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_a")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  # color rows
  opt_row_striping() %>%
  # width
  tab_options(table.width = pct(130))

Code

# custom function to summarize each column differently
compute_table_1_b <- function(x) {
  # if the column is numeric
  if (all(is.numeric(x))) {
    # and the sum is below 100
    if (sum(x) < 100) {
      # it is time
      # so we "circular" weighted average them
      y <- psych::circadian.mean(
        c(
          rep(x[1], 237),
          rep(x[2], 187)
        ),
        na.rm = T
      )
      y <- round(y, 1)
    } else {
      # otherwise it is just number
      # so we "just" weighted average them
      y <- weighted.mean(x, c(237, 187))
      y <- round(y, 1)
    }
  } else {
    y <- "\u00b1"
  }
  return(y)
}

Code

# shark and orca area = NA if no location data
data_dive <- data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
table_1_b <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      first(as.numeric(
        format(date_tz[DiveTypeName == "Benthic"],
          format = "%H"
        )
      ))
    ),
    hour_day_arrival = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      last(as.numeric(
        format(date_tz[DiveTypeName == "Benthic"],
          format = "%H"
        )
      ))
    ),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # tab_spanner(
  #   label = "per individual",
  #   columns = c(
  #     "n_dives_mean",
  #     "n_dives_plus_minus",
  #     "n_dives_sd",
  #     "n_benthic_dives_mean",
  #     "n_benthic_dives_plus_minus",
  #     "n_benthic_dives_sd",
  #     "hour_day_departure_mean",
  #     "hour_day_departure_plus_minus",
  #     "hour_day_departure_sd",
  #     "hour_day_arrival_mean",
  #     "hour_day_arrival_plus_minus",
  #     "hour_day_arrival_sd"
  #   )
  # )  %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_b")
  ) %>%
  # color rows
  opt_row_striping()

Code

# display first table
table_1_a

Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

	# seals	All dives			Benthic dives		% Benthic dives per area
		Total	w/ loc	w/o loc	Total	Proportion	Shark	Orca	Total	Overlap
Post Breeding	237	1,068,746	1,042,451	26,295	36,443	3.4%	79.8%	100.0%	100.0%	79.8%
Post Molting	187	2,207,361	2,032,528	174,833	44,594	2.0%	90.5%	100.0%	100.0%	90.5%
OVERALL	424	3,276,107	3,074,979	201,128	81,037	2.5%	87%	100%	100%	87%

Code

# display second table
table_1_b

Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

	# dives			# benthic dives			Hour departure			Hour arrival
	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
Post Breeding	4509.5	±	904.4	155.1	±	275.7	18.7	±	1.8	14.8	±	2.5
Post Molting	11804.1	±	1379.4	238.5	±	471.3	18.9	±	1.7	12.9	±	2.0
OVERALL	7726.7	±	1113.8	191.9	±	362	18.8	±	1.8	13.9	±	2.3

Notes

• Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

• Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

4.2 Figure 1

Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.

Code

# to choose an individual ~180 dives for 36 hours
data_dive %>%
  # rename date as datetime
  mutate(datetime = date) %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(datetime), datetime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(datetime), datetime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  )) %>%
  group_by(id, trip) %>%
  summarise(nb_dives = length(DiveNumber)) %>%
  arrange(nb_dives) %>%
  View()

Code

# import the csv export
tdr_data <- lapply(
  list.files(
    path = "../data/",
    pattern = "*_iknos_raw_data.csv",
    full.names = T
  ),
  function(x) {
    # extract id
    id_name <- last(strsplit(strsplit(x, split = "_")[[1]][1],
      split = "/"
    )[[1]])
    # load file
    dt <- readr::read_csv(x)
    # add id
    dt %>%
      mutate(id = as.numeric(id_name))
  }
) %>% bind_rows()

# choose the lag
nb_days <- 1.5

# compute the data
data_plot <- tdr_data %>%
  # convert time into date, a Posix
  mutate(
    datetime = matlab_to_posix(time),
    date = as.Date(datetime)
  ) %>%
  # add beginning and end date of the trip
  merge(
    .,
    data_dive %>%
      group_by(id, trip) %>%
      summarise(
        date_low = min(date),
        date_high = max(date),
        .groups = "drop"
      ),
    by = "id",
    all.x = TRUE
  ) %>%
  # identify what to keep
  mutate(to_keep = if_else(between(datetime, date_low, date_high), T, F)) %>%
  # keep only between the beginning and the end
  filter(to_keep) %>%
  # sort to make sure the difftime are being done correctly
  arrange(id, time) %>%
  # merge to get lat and lon
  merge(
    .,
    data_dive %>%
      group_by(id, date = as.Date(date)) %>%
      summarise(
        lat = median(Lat, na.rm = T),
        lon = median(Long, na.rm = T),
        .groups = "drop"
      ),
    by = c("id", "date"),
    all.x = TRUE
  ) %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(datetime), datetime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(datetime), datetime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  ))

# add bathy
setDT(data_plot)
setDT(data_dive)
data_dive[, datetime := date]
setkeyv(data_plot, c("id", "datetime"))
setkeyv(data_dive, c("id", "datetime"))
data_plot <- copy(data_dive) %>%
  .[, .(id, datetime, bathy, DiveTypeName)] %>%
  .[data_plot, roll = T]

# get sunrise sunset information
sun_data <- data_plot %>%
  select(date, lat, lon, id, phase) %>%
  unique() %>%
  ungroup() %>%
  mutate(getSunlightTimes(data = ., keep = c("sunrise", "sunset"))) %>%
  # fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1
  mutate(sunset = if_else( # if sunset is the next day
    abs(as.numeric(
      difftime(as.POSIXct(date), sunset, units = "days")
    )) > 1,
    # subtract 1
    sunset - (1 * 60 * 60 * 24),
    # otherwise leave it that way
    sunset
  ))

# trim if night beging before the trip or end after they arrive
sun_data <-
  merge(
    sun_data %>% select(date, lat, lon, id, phase, sunrise, sunset),
    data_plot %>%
      group_by(id, phase, trip) %>%
      summarise(
        date_low = min(datetime),
        date_high = max(datetime),
        .groups = "drop"
      ),
    by = c("id", "phase"),
    all.x = TRUE,
  ) %>%
  mutate(
    # make sure sunset is not outside the 36h period
    correct_sunset = if_else(
      # if not
      between(sunset, date_low, date_high),
      sunset,
      # if so, replace by one extremum
      if_else(sunset < date_low, date_low, date_high)
    ),
    # make sure sunrise is not outside the 36h period
    correct_sunrise = if_else(
      # if not
      between(sunrise, date_low, date_high),
      sunrise,
      # if so, replace by one extremum
      if_else(sunrise < date_low, date_low, date_high)
    )
  ) %>%
  # shift everything to start 2000-01-01 00:00:00
  mutate(
    correct_sunset_shift = correct_sunset - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    ),
    correct_sunrise_shift = correct_sunrise - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    )
  )

# plot
departure_dp <- data_plot %>%
  filter(phase == "Departure") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(datetime = datetime - (
    first(datetime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = datetime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Departure"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = datetime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("+6h", "+12h", "+18h", "+24h", "+30h"),
    position = "top"
  ) +
  scale_y_reverse() +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    strip.text.y = element_blank(),
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      # arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_line(
      color = "black",
      arrow = arrow(
        length = unit(0.2, "lines"),
        type = "closed",
        ends = "first"
      )
    ),
    axis.title.x = element_blank()
  )

arrival_dp <- data_plot %>%
  filter(phase == "Arrival") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(datetime = datetime - (
    first(datetime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = datetime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Arrival"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = datetime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("-30h", "-24h", "-18h", "-12h", "-6h"),
    position = "top"
  ) +
  scale_y_reverse() +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

fig_label_x <-
  ggplot(data.frame(l = "Time", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

Code

guide_area() / fig_label_x / (departure_dp | arrival_dp) +
  plot_layout(
    guides = "collect",
    heights = c(0.5, 0.7, 6)
  ) &
  theme(
    legend.position = "top",
    plot.margin = unit(c(0, 0, 0, 0), "npc")
  )

Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.

Notes

• TODO: Add three dots between Departure and Arrival

4.3 Figure 2

Let’s create a table specific for this figure that only contains for each individual:

• the first dive
• the last dive
• all benthic dives

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # then by individual
  group_by(id) %>%
  # keep only the first, last date, but also all benthic dives
  filter(date == min(date) |
    date == max(date) |
    DiveTypeName == "Benthic") %>%
  # ungroup
  ungroup() %>%
  # change colnames to make getSunlightTimes works
  mutate(
    date = as.Date(date_tz),
    lat = Lat,
    lon = Long
  ) %>%
  # calculate hour of sunset and hour of sunrise
  mutate(
    getSunlightTimes(
      data = .,
      keep = c("sunrise", "sunset"),
      tz = tz_lookup_coords(
        lat = Lat,
        lon = Long,
        method = "accurate"
      )
    )
  ) %>%
  # extract hours of sunset an sunrise
  mutate(
    h_sunrise = hour(sunrise),
    h_sunset = hour(sunset)
  )

Code

hour_time_2_radian <- function(x, high = 24, low = 0) {
  # https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles
  ptp <- high - low
  rad <- x * (2 * pi) / ptp
  rad <- (rad + pi) %% (2 * pi) - pi
  return(rad)
}

Code

# departure
test_departure_pm <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(DiveNumber == 1 & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)
test_departure_pb <- CircStats::r.test(
  data_windrose %>%
    group_by(id) %>%
    filter(DiveNumber == 1 & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)

# arrival
test_arrival_pm <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>% pull(time),
    units = "hours"
  ),
  units = "radians"
))
test_arrival_pb <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>% pull(time),
    units = "hours"
  ),
  units = "radians"
))

# # benthic
# test_benthic = r.test(conversion.circular(
#   circular(
#     data_windrose %>%
#       group_by(id) %>%
#       filter(DiveTypeName == "Benthic") %>%
#       mutate(time = hour(date_tz) +
#                minute(date_tz)/60 +
#                second(date_tz)/(60*60)) %>%
#       pull(time),
#     units = "hours"
#   ),
#   units = "radians"
# ))

# convert back to hours
conversion.circular(
  circular(test_departure_pm$r.bar,
    units = "radians"
  ),
  units = "hours"
)

Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
[1] 0.8608077

Code

conversion.circular(
  circular(test_departure_pb$r.bar,
    units = "radians"
  ),
  units = "hours"
)

Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
[1] 0.802508

Code

conversion.circular(
  circular(test_arrival_pm$r.bar,
    units = "radians"
  ),
  units = "hours"
)

Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
[1] 0.4919185

Code

conversion.circular(
  circular(test_arrival_pb$r.bar,
    units = "radians"
  ),
  units = "hours"
)

Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
[1] 0.1480673

Code

# conversion.circular(circular(test_benthic$r.bar,
#                              units = "rad"),
#                     units = "hours")

Code

# interesting circular time series
# https://campus.datacamp.com/courses/fraud-detection-in-r/introduction-motivation?ex=4

# for departure
windrose_departure_pm <- data_windrose %>%
  group_by(id) %>%
  filter(DiveNumber == 1 & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10, 15),
    facet_wrap = "origin",
    hour_sunrise = data_windrose %>%
      filter(DiveNumber == 1 &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      filter(DiveNumber == 1 &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  # geom_vline(xintercept = 6,
  #            linetype = "dashed") +
  geom_text(
    data = data.frame(
      x = rep(8, 3),
      y = c(6, 11, 16),
      text = c(5, 10, 15)
    ),
    aes(x = x, y = y, label = text),
    size = 2
  )

windrose_departure_pb <- data_windrose %>%
  group_by(id) %>%
  filter(DiveNumber == 1 & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10, 15),
    hour_sunrise = data_windrose %>%
      filter(DiveNumber == 1 &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      filter(DiveNumber == 1 &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
# geom_vline(xintercept = 6,
#            linetype = "dashed") +
# geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

# for arrival
windrose_arrival_pm <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10, 15),
    legend_position = "top",
    facet_grid = T,
    hour_sunrise = data_windrose %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

windrose_arrival_pb <- data_windrose %>%
  group_by(id) %>%
  filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n(), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10, 15),
    legend_position = "top",
    facet_wrap = "trip",
    hour_sunrise = data_windrose %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      median(.),
    hour_sunset = data_windrose %>%
      filter(date_tz == max(date_tz) &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      median(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

label_x <- ggplot(data.frame(l = "Time (hours)", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

# combine all
(windrose_departure_pm + windrose_arrival_pm) /
  (windrose_departure_pb + windrose_arrival_pb) /
  label_x +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  plot_layout(
    heights = c(2, 2, 0),
    guides = "collect"
  ) &
  # plot_annotation(tag_levels = c("A")) &
  # theme(plot.tag.position = c(-0.03, 0.8),
  #       plot.tag = element_text(size = 8, hjust = 0, vjust = 0)) +
  theme(legend.position = "top")

Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.

Note

• TODO: right now the shaded area is arbitrary defined, by I’ll compute the average sunset and sunrise time to be more accurate

Code

# departure PM
psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(DiveNumber == 1 & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

[1] 19.2026

Code

# departure PB
psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(DiveNumber == 1 & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

[1] 19.38994

Code

# arrival PM
psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

[1] 13.75121

Code

# arrival PB
psych::circadian.mean(
  data_windrose %>%
    group_by(id) %>%
    filter(date_tz == max(date_tz) & trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

[1] 15.68945

4.4 Figure 3

Code

# with or without bathymetry
with_bathy <- T

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_in.rds")) {
    trip_zoom_in <- readRDS("../export/background_ggoceanmap_zoom_in.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_in <- basemap(
      # limits = c(170, -110, 30, 59),
      limits = c(-155, -110, 30, 60),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_in)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_in)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_in <- trip_zoom_in +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_in$layers <- trip_zoom_in$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_in, "../export/background_ggoceanmap_zoom_in.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_in <- basemap(
    # limits = c(170, -110, 20, 59),
    limits = c(-155, -110, 30, 60),
    land.col = "black",
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_in)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_in)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_in <- trip_zoom_in +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_in$layers <- trip_zoom_in$layers[c(2, 1)]
}

Code

# get predator distribution *.kmz (extraction using google earth)
# for calculation
Shark_area_coord <- read_sf("../export/shark.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Orca_area_coord <- read_sf("../export/orca.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)

# for plots
Overlap_area_coord_to_display <- st_intersection(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Shark_area_coord_to_display <- st_difference(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)
Orca_area_coord_to_display <- st_difference(
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  select(lon = X, lat = Y)

Code

# add inside/outside each area
data_dive <- data_dive %>%
  mutate(
    shark_area = point.in.polygon(
      point.x = Long,
      point.y = Lat,
      pol.x = Shark_area_coord$lon,
      pol.y = Shark_area_coord$lat
    ),
    orca_area = point.in.polygon(
      point.x = Long,
      point.y = Lat,
      pol.x = Orca_area_coord$lon,
      pol.y = Orca_area_coord$lat
    ),
    overlap_area = as.numeric(if_else(shark_area == 1 &
      orca_area == 1, T, F)),
    total_area = as.numeric(if_else(shark_area == 0 &
      orca_area == 0, F, T))
  )

Code

# data_dive_longer <- pivot_longer(
#   data_dive %>%
#     # keep data with locations
#     filter(!is.na(Lat)),
#   c(
#     "shark_area",
#     "orca_area",
#     "overlap_area",
#     "total_area"
#   ),
#   names_to = "areas",
#   values_to = "in_out"
# )
#
# data_descriptive_ind <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, id) %>%
#   summarise(
#     n_tot = n(),
#     n_inside = length(id[in_out == 1]),
#     n_outside = length(id[in_out == 0]),
#     .groups = "drop"
#   ) %>%
#   mutate(
#     perc_inside = n_inside / n_tot,
#     perc_outside = n_outside / n_tot
#   )
#
# data_descriptive_area <- data_descriptive_ind %>%
#   group_by(areas) %>%
#   summarise(
#     mean_perc_inside = wtd.mean(perc_inside, n_tot),
#     mean_perc_outside = wtd.mean(perc_outside, n_tot),
#     sd_perc_inside = sqrt(wtd.var(perc_inside, n_tot)),
#     sd_perc_outside = sqrt(wtd.var(perc_outside, n_tot)),
#     # based from here: https://stats.stackexchange.com/questions/25895/computing-standard-error-in-weighted-mean-estimation
#     se_perc_inside = sqrt(wtd.var(perc_inside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#     se_perc_outside = sqrt(wtd.var(perc_outside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#   ) %>%
#   pivot_longer(
#     cols = ends_with("inside") |
#       ends_with("outside"),
#     names_to = c("stat", "in_out"),
#     names_pattern = "(.*)_perc_(.*)"
#   ) %>%
#   pivot_wider(
#     id_cols = c("areas", "in_out"),
#     names_from = "stat",
#     values_from = "value"
#   ) %>%
#   mutate(areas = factor(areas,
#                         level = c(
#                           "overlap_area",
#                           "orca_area",
#                           "shark_area",
#                           "total_area"
#                         )
#   ))
#
# data_stat_area <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, in_out) %>%
#   summarise(nb_area = n(), .groups = "drop_last") %>%
#   mutate(nb_total = sum(nb_area)) %>%
#   summarise(rstatix::prop_test(x = nb_area, n = nb_total))

Code

# inside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside", -mean, mean),
#     sd = if_else(in_out == "outside", -sd, sd),
#     se = if_else(in_out == "outside", -se, se)
#   ) %>%
#   filter(in_out == "inside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(
#     show.legend = F,
#     mapping = aes(
#       x = mean,
#       y = areas,
#       fill = in_out),
#     alpha = 0.3
#   ) +
#   scale_fill_manual(values = c("grey70")) +
#   geom_errorbar(mapping = aes(xmin = -se,
#                               xmax = se),
#                 # reduce size of horizontal bar of the errorbar
#                 width = 0.05) +
#   coord_flip(xlim = c(0, max(data_descriptive_area$mean))) +
#   scale_x_continuous(
#     # breaks = c(-0.6,-0.4,-0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0,0.1), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x)
#         (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   theme(
#     panel.spacing.y = unit(0, "mm"),
#     panel.grid.major = element_blank(),
#     panel.grid.minor = element_blank(),
#     panel.background = element_blank(),
#     axis.line.y = element_line(arrow = arrow(
#       angle = 20,
#       length = unit(.10, "inches"),
#       type = "closed"
#     )),
#     axis.title = element_blank(),
#     axis.ticks.length.x = unit(0, units = "mm"),
#     axis.text.x = element_blank()
#   )
#
# outside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside",-mean, mean),
#     sd = if_else(in_out == "outside",-sd, sd),
#     se = if_else(in_out == "outside",-se, se)
#   ) %>%
#   filter(in_out == "outside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(show.legend = F, mapping = aes(x = mean,
#                                           y = areas,
#                                           fill = in_out),
#                                           alpha = 0.3) +
#   scale_fill_manual(values = c("grey30")) +
#   geom_errorbar(
#     mapping = aes(
#       xmin = -se,
#       xmax = se
#     ),
#     # reduce size of horizontal bar of the errorbar
#     width = 0.05
#   ) +
#   coord_flip(xlim = c(-max(data_descriptive_area$mean), 0)) +
#   scale_x_continuous(
#     # breaks = c(-0.6, -0.4, -0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0.1,0), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   labs(y = "Areas") +
#   theme(panel.spacing.y = unit(0, "mm"),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.background = element_blank(),
#         axis.line.y = element_line(arrow = arrow(
#           angle = 20,
#           length = unit(.10, "inches"),
#           type = "closed",
#           ends = "first"
#         )),
#         axis.ticks.length.x = unit(0,units = "mm"),
#         axis.title.y = element_blank())
#
# label_x_bar_plot = ggplot(data.frame(l = "Proportion of benthic dives",
#                                      x = 1,
#                                      y = 1)) +
#   geom_text(aes(x, y, label = l, angle = 90)) +
#   theme(aspect.ratio = 10,
#         axis.title = element_blank(),
#         axis.text = element_blank(),
#         axis.ticks = element_blank()) +
#   theme_void()

Code

# color palette
# https://davidmathlogic.com/colorblind/#%23D81B60-%231E88E5-%23FFC107-%23004D40
# combine areas
areas_coord <- rbind(
  Orca_area_coord_to_display %>% mutate(Predator = "orca"),
  Shark_area_coord_to_display %>% mutate(Predator = "shark"),
  Overlap_area_coord_to_display %>% mutate(Predator = "overlap")
) %>%
  mutate(Predator = factor(Predator, levels = c("orca", "shark", "overlap")))

# data use to compute kernel density estimation
df_kernel_dens_benthic <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

# transform data into sf object
df_kernel_dens_sf_benthic <- st_as_sf(
  df_kernel_dens_benthic,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf_benthic <- group_by(df_kernel_dens_sf_benthic, origin)

# kernel density estimation
df_kernel_dens_sf_kde_benthic <- st_kde(df_kernel_dens_sf_benthic,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf_benthic)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf_benthic)[, 2]
    )
  ) / 4)^2
)

# https://loading.io/color/feature/BuGn-9/
cols <- c(
  "99%" = "#ccece6",
  "95%" = "#99d8c9",
  "80%" = "#41ae76",
  "50%" = "#00441b"
)

p1 <- trip_zoom_in +
  # remove legend
  guides(fill = "none") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # legend
  labs(fill = "Probs") +
  # legend at the top
  theme(legend.position = "top") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# Extract the size legend - leg1
leg1 <- gtable_filter(ggplot_gtable(ggplot_build(p1)), "guide-box")

p2 <- trip_zoom_in +
  # tweak legend
  guides(fill = guide_legend(
    override.aes = list(
      colour = "grey60"
    )
  )) +
  theme(
    legend.key = element_rect(color = "grey60"),
    legend.text = element_text(colour = "grey20"),
    legend.title = element_text(colour = "grey20"),
    legend.justification = c("right", "top"),
    legend.position = c(0.99, 0.99),
    legend.background = element_rect(fill = alpha("white", 0))
  ) +
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # remove legend
  guides(fill = "none") +
  # # no display alpha
  # guides(alpha = "none") +
  # new scale
  new_scale_fill() +
  # add spatial polygon
  geom_spatial_polygon(
    data = areas_coord,
    aes(
      x = lon,
      y = lat,
      fill = Predator
    ),
    alpha = 0.35,
    crs = 4326
  ) +
  # color palette
  scale_fill_manual(values = c("darkred", "#E1E826", "#E66100")) +
  scale_color_manual(values = c("darkred", "#E1E826", "#E66100")) +
  # remove legend for areas
  guides(fill = "none", col = "none") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# reordering layers (continent on top)
p2$layers <- if (with_bathy) {
  p2$layers[c(1, 2, 5, 4, 3)]
} else {
  p2$layers[c(1, 4, 3, 2)]
}

# convert into ggplot
predator_map <- as_ggplot(leg1) / p2 +
  plot_layout(
    ncol = 1,
    heights = c(1, 12)
  )

Code

predator_map

Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.

Code

# # layout
# layout <- "
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# #AAAAAAAAAAAAAAAAAAAAAAA
# BCCCCCCCCCCCCCCCCCCCCCCC
# BDDDDDDDDDDDDDDDDDDDDDDD
# "
#
# # predator_map$layers[[length(predator_map$layers)]]$aes_params$fill_new_new = "white" #cdaa6d
# # predator_map$layers[[length(predator_map$layers)]]$aes_params$colour = "black"
#
# # display
# (predator_map + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (label_x_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (inside_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   (outside_bar_plot + theme(plot.margin = unit(c(0, 0, 0, 0), "pt"))) +
#   plot_layout(design = layout) #+
#   #plot_annotation(tag_levels = list(c("A","","B","")))

Notes

• ?

4.5 Figure 4

Code

# main plot
bathy_prop <- data_dive %>%
  # keep only negative bathy
  filter(bathy < 0) %>%
  # and remove outliers
  filter(bathy > -6000) %>%
  # positive bathy
  mutate(bathy = abs(bathy)) %>%
  # create class of bathymetry
  mutate(bathy_class = cut(
    bathy,
    seq(0, 6000, 400),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(bathy_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per bathy_class
  group_by(bathy_class) %>%
  # to get the Proportion of different dive types per bathy_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = bathy_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Bathymetry  (m)",
    y = "Dive type proportion (%)"
  ) +
  scale_fill_viridis(
    "Dive Type",
    option = "plasma",
    discrete = T,
    direction = -1
  ) +
  # scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))
  # position legend at the top
  theme(
    legend.position = "top",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  )

# # second plot
# bathy_hist <- data_dive %>%
#   # keep only negative bathy
#   filter(bathy < 0) %>%
#   # and remove outliers
#   filter(bathy > -6000) %>%
#   # create class of bathymetry
#   mutate(bathy_class = fct_rev(cut(
#     bathy,
#     seq(-6000, 0, 400),
#     ordered_result = T,
#     dig.lab = 4
#   ))) %>%
#   # calculate by bath_class and animal
#   group_by(bathy_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = bathy_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))

Code

# main plot
dist_coast_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_coast > 0) %>%
  # remove outliers
  filter(dist_coast * 1000 * 111 <= 1900) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    # convert decimal degree/1000
    dist_coast * 1000 * 111,
    seq(0, 1900, 100),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from the coast (km)",
    y = "Dive type proportion"
  ) +
  scale_fill_viridis(
    "Dive Type",
    option = "plasma",
    discrete = T,
    direction = -1
  ) +
  # scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  )

# # the second plot
# dist_coast_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_coast * 1000 * 111 <= 1900) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     # convert decimal degree/1000
#     dist_coast * 1000 * 111,
#     seq(0, 2000, 100),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# the actual plot
# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))

Code

# the main plot
dist_dep_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_dep > 0) %>%
  # remove outliers
  filter(dist_dep / 1000 < 4200) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    dist_dep / 1000,
    seq(0, 4200, 300),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from departure (km)",
    y = "Dive type proportion"
  ) +
  scale_fill_viridis(
    "Dive Type",
    option = "plasma",
    discrete = T,
    direction = -1
  ) +
  # scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  )

# # the second plot
# dist_dep_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_dep / 1000 < 4200) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     dist_dep / 1000,
#     seq(0, 4200, 300),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))

Code

# display
ylab <- bathy_prop$labels$y
bathy_prop$labels$y <- dist_coast_prop$labels$y <- dist_dep_prop$labels$y <- " "
(guide_area() + bathy_prop + dist_coast_prop + dist_dep_prop) +
  plot_annotation(tag_levels = "A") +
  plot_layout(nrow = 5, guides = "collect")
grid::grid.draw(grid::textGrob(ylab, x = 0.06, y = 0.55, rot = 90))

Figure 4: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).

Notes

• TODO: change the color code to match the color we attributed to benthic dives (greens?)

4.6 Figure 5

Code

# the last plot
prop_at_sea <- data_dive %>%
  group_by(id) %>%
  # time difference between the first date and the others (telling R to take diff b/w first date and all other dates)
  mutate(nb_days_departure = trunc(as.numeric(difftime(
    date, first(date),
    units = "days"
  )))) %>%
  # group by day/dive_type/seal
  group_by(nb_days_departure, DiveTypeName, id) %>%
  # count number of dives
  summarise(nb_daily_dives_type = n(), .groups = "drop") %>%
  # group by day/seal
  group_by(nb_days_departure, id) %>%
  # count the total number of dives per day/seal
  mutate(nb_daily_dives = sum(nb_daily_dives_type)) %>%
  # order by seals/date
  arrange(id, nb_days_departure) %>%
  # perform our calculation per seal
  group_by(id) %>%
  # calculate of the Proportion of time since departure
  mutate(perc_time_at_sea = round(
    nb_days_departure / max(nb_days_departure) * 100, 1
  )) %>%
  # filter on benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # add class of Proportion of time at sea
  mutate(day_class = cut(perc_time_at_sea, seq(0, 100, 5),
    include.lowest = T
  )) %>%
  # by class of Proportion of time at sea
  group_by(day_class) %>%
  # calculate the proportion of daily benthic dives
  summarise(
    nb_benthic_class = sum(nb_daily_dives_type),
    nb_total_class = sum(nb_daily_dives),
    prop_class = nb_benthic_class / nb_total_class
  ) %>%
  # plot
  ggplot(
    .,
    aes(x = day_class, y = prop_class)
  ) +
  geom_bar(stat = "identity", fill = "#008c49") +
  guides(x = guide_axis(angle = 45)) +
  scale_y_continuous(
    limits = c(0, 1),
    # breaks=seq(0,0.6,0.1),
    labels = function(x) {
      scales::percent(abs(x), 1)
    }
  ) +
  annotate("text", x = 5.3, y = 0.4, label = "Departure") +
  annotate("text", x = 15, y = 0.5, label = "Arrival") +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(4, 17),
  #     x2 = c(1.5, 19.5),
  #     y1 = c(0.4, 0.5),
  #     y2 = c(0.31, 0.55)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = 0.3
  # ) +
  geom_curve(
    data = tibble(
      x1 = c(4),
      x2 = c(1),
      y1 = c(0.4),
      y2 = c(0.33)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = 0.3
  ) +
  geom_curve(
    data = tibble(
      x1 = c(16),
      x2 = c(19.7),
      y1 = c(0.5),
      y2 = c(0.57)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = -0.56
  ) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  ) +
  labs(
    x = "Time spent at sea (%)",
    y = "Proportion of benthic dives"
  )

Code

# display
prop_at_sea

Figure 5: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.

Notes

• ?

4.7 Supplementary

4.7.1 Table 1

Code

# shark and orca area = NA if no location data
data_dive <- data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
data_dive %>%
  group_by(year, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" & !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      first(as.numeric(format(date_tz[DiveTypeName == "Benthic"],
        format = "%H"
      )))
    ),
    hour_day_arrival = ifelse(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      last(as.numeric(format(date_tz[DiveTypeName == "Benthic"],
        format = "%H"
      )))
    ),
    .groups = "drop"
  ) %>%
  group_by(year) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt() %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  tab_spanner(
    label = "per individual",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # rename columns
  cols_label(
    year = "Year",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark area",
    n_benthic_dives_in_orca_proportion = "Orca area",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # color rows
  opt_row_striping() %>%
  tab_options(table.width = pct(180))

											per individual
Year	# seals	All dives			Benthic dives		% Benthic dives				# dives			# benthic dives			Hour departure			Hour arrival
		Total	w/ loc	w/o loc	Total	Proportion	Shark area	Orca area	Total	Overlap	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
2004	26	280,155	215,452	64,703	9,291	3.3%	79.0%	100.0%	100.0%	79.0%	10775.2	±	3031.8	357.3	±	575.3	17.7	±	1.5	14.4	±	1.4
2005	35	281,636	226,566	55,070	16,413	5.8%	57.8%	100.0%	100.0%	57.8%	8046.7	±	3704.9	468.9	±	958.3	18.7	±	1.4	19.2	±	2.0
2006	29	219,773	181,667	38,106	5,219	2.4%	72.5%	100.0%	100.0%	72.5%	7578.4	±	3613.8	180.0	±	232.1	18.7	±	1.5	11.3	±	1.6
2007	30	235,691	235,691	0	5,829	2.5%	91.4%	100.0%	100.0%	91.4%	7856.4	±	3880.7	194.3	±	366.8	18.5	±	1.6	15.4	±	2.1
2008	35	252,268	239,800	12,468	3,460	1.4%	99.9%	100.0%	100.0%	99.9%	7207.7	±	3557.8	98.9	±	135.7	13.0	±	2.2	12.6	±	2.0
2009	24	149,429	135,010	14,419	3,646	2.4%	94.4%	100.0%	100.0%	94.4%	6226.2	±	3076.3	151.9	±	188.2	15.3	±	1.8	18.1	±	1.8
2010	35	248,956	248,956	0	3,702	1.5%	93.0%	100.0%	100.0%	93.0%	7113.0	±	3475.9	105.8	±	81.0	20.6	±	1.6	11.9	±	1.8
2011	32	244,228	232,301	11,927	3,883	1.6%	85.3%	100.0%	100.0%	85.3%	7632.1	±	3842.4	121.3	±	151.4	22.7	±	1.7	7.8	±	1.7
2012	33	256,295	256,295	0	4,233	1.7%	82.2%	100.0%	100.0%	82.2%	7766.5	±	4090.7	128.3	±	152.0	20.1	±	1.3	13.4	±	2.6
2013	31	265,526	265,526	0	5,149	1.9%	87.0%	100.0%	100.0%	87.0%	8565.4	±	4318.0	166.1	±	218.6	15.2	±	2.2	6.9	±	1.7
2014	22	141,911	141,911	0	6,679	4.7%	90.5%	100.0%	100.0%	90.5%	6450.5	±	3424.0	303.6	±	473.9	17.0	±	1.7	17.7	±	1.7
2015	30	226,090	226,090	0	4,380	1.9%	99.3%	100.0%	100.0%	99.3%	7536.3	±	3956.1	146.0	±	91.3	17.6	±	1.1	18.9	±	2.0
2016	27	215,401	215,401	0	4,003	1.9%	93.0%	100.0%	100.0%	93.0%	7977.8	±	4315.7	154.0	±	101.5	22.1	±	1.6	13.1	±	1.7
2017	10	75,326	75,326	0	1,648	2.2%	65.8%	100.0%	100.0%	65.8%	7532.6	±	3664.9	164.8	±	162.3	5.4	±	1.6	0.1	±	1.0
2018	18	150,599	146,164	4,435	2,733	1.8%	92.0%	100.0%	100.0%	92.0%	8366.6	±	3650.2	160.8	±	89.8	11.2	±	2.2	10.4	±	2.2
2019	7	32,823	32,823	0	769	2.3%	97.8%	100.0%	100.0%	97.8%	4689.0	±	1342.1	109.9	±	65.8	15.7	±	1.6	20.1	±	2.6

Notes

• lmk!

4.7.2 Figure 1

Let’s first download the map

Code

# with or without bathymetry
with_bathy <- F

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_out.rds")) {
    trip_zoom_out <- readRDS("../export/background_ggoceanmap_zoom_out.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_out <- basemap(
      limits = c(170, -110, 30, 59),
      # limits = c(-155, -110, 30, 60),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_out)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_out)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_out <- trip_zoom_out +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_out$layers <- trip_zoom_out$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_out, "../export/background_ggoceanmap_zoom_out.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_out <- basemap(
    limits = c(170, -110, 20, 59),
    # limits = c(-155, -110, 30, 60),
    # land.col = 'black',
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_out)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_out)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_out <- trip_zoom_out +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_out$layers <- trip_zoom_out$layers[c(2, 1)]
}

Tip

This part is mostly based on https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html. The author of the package confirms that his kernel density estimation is calculated using “meaningful” units such as UTM ¹, since the input is in a simple feature format.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat)) %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf <- group_by(df_kernel_dens_sf, origin)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
to_plot <- trip_zoom_out +
  # rename axis
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # facet by origin
  facet_wrap(. ~ origin)

# reorder layers
to_plot$layers <- if (with_bathy) {
  to_plot$layers[c(1, 2, 4, 3)]
} else {
  to_plot$layers[c(1, 3, 2)]
}

# display
to_plot

Figure 6: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probs represents different levels of probability density.

Notes

• lmk!

4.8 Extra

Code

# import the csv export
data_plot <- readr::read_csv("../data/benthic_dive_extract.csv")

# plot
data_plot %>%
  ggplot(aes(x = time, y = CorrectedDepth)) +
  geom_path() +
  scale_y_reverse() +
  theme_minimal() +
  theme(
    axis.line.y = element_line(arrow = arrow(
      type = "closed",
      ends = "first",
      length = unit(0.1, "inches")
    )),
    axis.line.x = element_line(arrow = arrow(
      type = "closed",
      ends = "last",
      length = unit(0.1, "inches")
    ))
  ) +
  labs(
    x = "Time",
    y = "Depth (m)"
  ) +
  # geom_text(data = data.frame(
  #          xpos = median(data_plot$time),
  #          ypos = 30,
  #          hjustvar = 0,
  #          vjustvar = 0,
  #          annotateText = paste0("Corner index: ", round(unique(data_plot$corner_index),2),
  #                       "\n",
  #                       "Benthix index: ", round(unique(data_plot$benthic_index)),2)),
  #           aes(x=xpos,y=ypos,hjust=hjustvar,
  #               vjust=vjustvar,label=annotateText)) +
  scale_x_datetime(position = "top") +
  labs(subtitle = "ID: 2011034 - Dive: 29")

Footnotes

1. • https://gis.stackexchange.com/questions/64638/how-to-fill-in-the-parameters-for-kernel-density-estimation
   • https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1020&context=siv
   • https://gis.stackexchange.com/questions/50692/choosing-projection-for-rasterization-of-random-latitude-and-longitude-data-in-n/50791#50791